Source Tracing of Leishmania donovani in Emerging Foci of Visceral Leishmaniasis, Western Nepal

We sequenced Leishmania donovani genomes in blood samples collected in emerging foci of visceral leishmaniasis in western Nepal. We detected lineages very different from the preelimination main parasite population, including a new lineage and a rare one previously reported in eastern Nepal. Our findings underscore the need for genomic surveillance.


Laboratory procedures DNA extraction and sequencing
200µl of blood mixed with 200µl DNA/RNA Shield was used to extract DNA with QIAamp DNA Mini Kit (Qiagen) according to manufacturer's instructions with the following modifications: 30µl Proteinase K and 300µl ethanol were used instead of 20µl and 200µl, respectively.DNA concentration was verified using the Qubit broad-range DNA quantification kit (Thermo Fisher Scientific), and the % of Leishmania DNA in the samples was estimated using qPCR as described previously (1).In our previous study (1), we demonstrated that a percentage of Leishmania DNA of 0.006% was found to be the lowest limit for suitable analysis of genome diversity.In present study, the Leishmania % in the selected samples was 0.13, 0.16 and 0.03 in samples 022, 023, and 024, respectively.SureSelect (Agilent Technologies) was used to capture Leishmania genomic DNA following standard SureSelect XTHS Target Enrichment system protocol for Illumina Multiplexed Sequencing platforms.Prior the genome capture, DNA was concentrated using AMPure XP beads (Beckman, Coulter) to obtain ≈10ng of total genomic DNA in 7 µl that was subjected to enzymatic fragmentation using SureSelect XT HS Fragmentation kit (Agilent Technologies).Custom designed oligonucleotide baits were used at 1:10 stock dilution.
Sequencing was conducted on the Illumina NovaSeq platform using 2x150 bp paired reads at GenomeScan (Netherlands), for which 51.49, 51.59 and 50.35 million raw reads were obtained for sample 022, 023 and 024 respectively (Appendix Table 2).

Bioinformatic procedures
In addition to the newly sequenced data, additional sequencing data were obtained from previous publications: i) genomes describing the population structure of L. donovani in Nepal, India and Bangladesh (2), ii) sequencing data obtained from samples from Nepal using the SureSelect technology (1), similar to the approach used for the three outbreak samples in this report, iii) three genomes originating from Sri Lanka (3,4), and iv) the L. infantum sequencing data submitted under the accession number ERR1913337 (5).All publicly available sequencing data were downloaded using the SRAtoolkit software.
The reads were mapped to the reference genome L. donovani available at NCBI (accession number GCF_000227135.1) using BWA (version 0.7.17 (6)) with a seed length set to 50 17 (6).Only properly paired reads with a mapping quality higher than 30 were selected using SAMtools (7).Duplicate reads were removed using the RemoveDuplicates command in the Picard software (version 2.22.4,http://broadinstitute.github.io/picard/).SNP calling was performed using the Genome Analysis ToolKit (GATK) (8) pipeline (version 4.1.4.1) following the GATK best practices approach: 1) GATK HaplotypeCaller enabling the production of GVCF formatted files, 2) GVCF files of all samples were combined using the GATK CombineGVCF command, 3) genotyping was performed via the GATK GenotypeGVCF command, and 4) filtering of the SNPs and indels was carried out following the "best practices" approach as suggested on the GATK support site using the SelectVariants and VariantFiltration commands.
Regions in the vcf-file corresponding to known drug resistance markers were selected using BCFtools, and visualized using the pheatmap function in R.
Phylogenetic trees were constructed using RAxML (9).First, the VCF files containing biallelic SNPs were selected using BCFtools (7) and were converted to Phylip format using the vcf2phylip.pyscript (https://github.com/edgardomortiz/vcf2phylip).RAxML was then executed with the GTR+G substitution model, using 1000 bootstrap replicates.The L. infantum JPCM5 or the L. donovani LV9 genome was employed as an outgroup.The resulting phylogenetic trees were visualized using ggtree (10) for rooted phylogenetic trees and SplitsTree (11) for unrooted phylogenetic networks.

Results from the analyses of specific genes reported to be involved in drug resistance
We selected 10 loci that were previously shown to be involved in L. donovani resistance to known antileishmanial drugs: For each of the three new L. donovani genomes, the sequence of the 10 loci was studied in detail.The 10 loci were well-covered and in 8 out of the 10 genes we found at least one homozygous single nucleotide polymorphism (SNP), which results in a missense mutation or a frameshift in at least 2 out of 3 samples from new emerging loci (Appendix Figure 2).No significant changes were observed for LdRos3 and SMT.

Functional differences previously reported between ISC1 and CG isolates
Noteworthy, all results here compiled concern the analysis of isolated and cultivated parasites.
In a first study, we demonstrated that CG parasites were intrinsically more tolerant to trivalent antimonials than ISC1 ones.This phenotype was driven by the amplification of a locus containing MRPA, a gene involved in Sb III sequestration (13) In a second study, we made an integrated genomic and metabolomic profiling of ISC1 vs CG isolates (17).We found several genomic differences including SNPs, CNV and small indels in genes coding for known virulence factors, immunogens and surface proteins.With respect to the metabolome, we found differences in several functional groups and pathways, essentially: (iii) Nucleotide salvage pathway.This pathway is essential, since Leishmania cannot synthesize the purine ring de novo and is therefore dependent on salvaging these from host purines.Our previous results suggested that ISC1 parasites might be better at salvaging nucleotides from their environment.
In a third study, we experimentally demonstrated that ISC1 and CG strains are developing similarly in natural ISC vector Phlebotomus argentipes, suggesting that P. argentipes is a fully competent vector for ISC1 parasites (18).
Altogether, these experimental studies demonstrate differences between ISC1 and CG in antimonial susceptibility and predict major functional differences, including virulence.Taking into account that ISC1 can easily be transmitted by P. argentipes, particular attention is required to monitor the fate of ISC1-related population in the region, especially in a post-VL elimination context.Genomic surveillance can clearly be done with the approach here used, but this should be complemented by phenotyping of the new detected L. donovani variants.

( i )
Lipid metabolism, with 19 glycerophospholipids (GPLs) showing significantly different levels between both groups: GPLs are involved in a wide array of cellular functions including host cell infection (ii) Urea cycle.In ISC1 versus CG we detected a higher concentration of citrulline and a lower concentration of argininosuccinate: Mutants for argininosuccinate synthase genes have shown a lower virulence than WT parasites.